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Exact ground states of three-dimensional random field Ising magnets (RFIM) with Gaussian dis- 
tribution of the disorder are calculated using graph-theoretical algorithms. Systems for different 
strengths h of the random fields and sizes up to iV = 96'^ are considered. By numerically differ- 
entiating the bond-energy with respect to ft a specific-heat like quantity is obtained, which does 
not appear to diverge at the critical point but rather exhibits a cusp. We also consider the effect 
of a small uniform magnetic field, which allows us to calculate the T — susceptibility. From a 
finite-size scaling analysis, we obtain the critical exponents p — 1.32(7), a = —0.63(7), rj — 0.50(3) 
and find that the critical strength of the random field is = 2.28(1). We discuss the significance 
of the result that a appears to be strongly negative. 



PACS numbers; 75.50.Lk, 05.70.Jk, 75.40. Mg, 77.80. Bh 



I. INTRODUCTION 

The random field Ising model Q has been extensively 
studied both because of its interest as a "simple" 

frustrated system and because of its relevance to experi- 
ments, especially those on the diluted antiferromagnet in 
a uniform field |p|. The RFIM Hamiltonian is given by 



H — — J ^ SiSj — ^ h^Si, 

{hi) i 



(1) 



where the S'i = ±1 are Ising spins, J is the interaction 
energy between nearest neighbors, and hi is the random 
field. The values hi are independently distributed accord- 
ing a Gaussian distribution with mean and standard 
deviation h, i.e. the probability distribution is 



P{h, 



1 



rhrh 



exp 



2K^ 



(2) 



We shall consider three-dimensional lattices with periodic 
boundary condition and N = spins. 

A sketch of the phase boundary is shown in Fig. |l|. At 
low values of the random field and temperature T, the 
system is in a ferromagnetic phase, and at high temper- 
atures or random fields, the system is paramagnetic. 

In this paper we shall be interested in the values of 
the critical exponents along the phase boundary. The 
random field is a relevant perturbation at the pure (i.e. 
h = 0) fixed point, and the random-field fixed point is 
at T = 1^,0. Hence, the critical behavior is the same 
everywhere along the phase boundary in Fig. ^ (assum- 
ing that the transition is always second order) except for 
h = 0. We can therefore determine the critical behavior 
by staying at T = and crossing the phase boundary at 
h = he, see Fig. |l|, which is convenient, because we can 
determine the ground states of large lattices exactly us- 
ing efficient optimization algorithms |]8|"|lT|, as discussed 
in Sec. ||. This has the advantage that one can study 
much larger systems than it is possible in Monte Carlo 
simulations, and, for each realization, there are no sta- 
tistical errors or equilibration problems. 




FIG. 1. A sketch of the phase boundary of the random 
field Ising model. The ferromagnetic phase is denoted by "F" 
and the paramagnetic phase by "P" . The critical value of the 
random field at T = is denoted by he- The lines with arrows 
at both ends indicate the path followed by varying J for some 
fixed value of h and T. 



Using these ground state techniques, most of the criti- 
cal exponents have been determined with some precision; 
for a thorough recent study see Ref . . Most of these 
exponents are consistent with scaling relations. How- 
ever, as we shall discuss in Sec. 0, those scaling relations 
predict a specific-heat exponent a close to zero, while 
Monte Carlo data on fairly small sizes (L < 16) find 
ajv = —0.45 ± 0.05, where v is the correlation length 
exponent (which has a value slightly greater than unity, 
as discussed in Sees. IV and 0). Interestingly, experi- 
ments find a logarithmic divergence, corresponding 
to a specific heat exponent a = 0, as expected from scal- 
ing. 



In order to try to resolve this puzzle, we calculate here 
the specific heat exponent for the RFIM using much 
larger sizes (L < 96) than in the Monte Carlo work 
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p3[ , by using optimization methods to determine exact 
ground states. We also find a strongly negative value 
for a, oijv = —0.48 ± 0.05, consistent with the earlier 
Monte Carlo data jl^, but in disagreement with expe- 
riment and apparently in violation of scaling. In Sec. ^ 
we will discuss possible ways round this discrepancy. In 
addition, we determine the susceptibility, which, to our 
knowledge, has not been directly computed before using 
ground-state methods. Our results are consistent with 
earlier calculations. 



The total energy per spin, i?, is given by 
E = JEj + hEh, 
where the "field energy" E^ is given by 



E, 



dF 
'dh 



1 



E 



{Si. 



(4) 



(5) 



II. NUMERICAL TECHNIQUES 

We used well known algorithms P-pl]] from graph the- 
ory |p"5|"p7t to calculate the ground state of a system 
at given random- field strength h. To implement them 
we applied some algorithms from the LEDA library . 
The calculation works by transforming the system into 
a network [ p^ , and calculating the maximum flow in 
polynomial time [^o|-^. The first results of applying 
these al gor ithms to random-field systems can be found 
in Ref. [2^]. In Ref. |26| these methods were applied 
to obtain the exponents for the magnetization, the dis- 
connected susceptibility and the correlation length from 
ground-state calculations up to size L = 80. Other ex- 
act ground-state calculation of the RFIM can be found 
in Refs. l?^- |2g| , [l^ . Note that in cases where the ground- 
state is degenerate it is possible to calculate all the 
ground-states in one sweep [0, see also Refs. 
For the RFIM with a Gaussian distribution of fields, the 
ground state is non-degenerate, except for a two-fold de- 
generacy at certain values of the randomness, where the 
ground state changes, see Sec. [II, so it is sufficient to 
calculate just one ground state. 



III. QUANTITIES OF INTEREST 

In zero random field, the specific-heat exponent is ob- 
tained from the singularity in the second derivative of the 
free energy with respect to temperature. More generally 
it is determined from the singularity obtained by varying 
a parameter which crosses the phase boundary from the 
paramagnetic phase to the ferromagnetic phase. From 
Fig. |l| we see that this can be conveniently accomplished 
by keeping the ratio of h/ J to T/ J fixed, i.e. by varying 
J. The first derivative of the free energy (per spin) F 
with respect to J, which we call the "bond energy" 
is given by 



Ej 



dF 



(3) 



where (• • •) is a thermal average, and the sum is over 
nearest-neighbor pairs. Ej has an energy-like singularity 
in the vicinity of the phase boundary. For /i = it is 
precisely the energy, apart from an overall factor of J. 
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FIG. 2. Bond energy per spin, Ej, defined in Eq. (^), for 
two L = 8 samples as a function of the random-field strength 
h. 

Having differentiated analytically with respect to J, 
we now set J = 1, consider T — Q only, and obtain a 
specific heat-like quantity by differentiating Ej numeri- 
cally with respect to the random field h. We emphasize 
that it is not necessary to vary the temperature in order 
to observe the specific heat singularity. To observe this 
singularity the direction in which the phase boundary is 
crossed must have a projection on to the correct scaling 
field, which means that the phase boundary should not be 
approached tangentially. The angle at which the phase 
boundary is approached will affect the size of corrections 
to scaling by mixing in a varying amount of irrelevant op- 
erators, but the asymptotic behavior will always be the 
same (as long as the approach is not tangential). 

To avoid confusion we point out that the role taken 
by the free energy at finite-T is played by the energy 
at T = 0, since the two are equal in this limit. More 
precisely, the energy singularity at T = has the form 
e^~", where e is the deviation from criticality, which is 
the same as the free- energy singularity at a finite-T tran- 
sition. At finite-T, the energy and entropy each have a 



2 



stronger singularity, of the form e^~°', but with opposite 
signs such that this singularity cancels in the free en- 
ergy, F = E — T S . A analogous cancellation occurs at 
T = 0, but between Ej and Eh since both Ej and Eh 
have singularities with exponent 1 — a but with ampli- 
tudes of opposite sign such that this singularity cancels 
in the total energy. To see this note that from Eq. 



dE _ jdEj ^ ^dEh 
dh dh dh 



Eh 



(6) 



However, at T = where F = E, we have dE/dh = Eh, 
and so, in this limit. 



J- 



dh 



(7) 



Hence, if Eh ^ \h ~ hc\^-", then dEh/dh and dEj/dh 
each have singularities of the form /ic|^", but with op- 
posite signs such that this singularity cancels in dE/dh. 
We have verified that this cancellation occurs in our nu- 
merical data. From Eqs. ^ and (^, dE/dh has the same 
singularity as Eh, i.e. \h — /icT so E ^ \h — hc\'^~°', 
as stated above. 

We use the same set of random fields for different val- 
ues of h and scale them all by a fixed overall factor. More 
precisely we take hi — eih, where the u are chosen from 
a Gaussian distribution with standard deviation unity, 
and are the same for all values of h. We use a first- 
order finite difference to determine the derivative of Ej 
numerically and, since this is a more accurate represen- 
tation of the derivative at the midpoint of the interval 
than at either endpoint, the "specific heat", C, at T = 
is defined to be 



C 



hi 



[Ej{hi)]h-[Ej{h2)]h 
hi — h2 



(8) 



where hi and /12 are two "close-by" values of h, and 
[■ ■ ■]h denotes an average over random-field configura- 
tions, which is carried out (approximately) by repeating 
the calculation for -/Vgamp independent realizations (sam- 
ples) of the random fields e^. We choose a sufficiently 
fine mesh of random-field values that the resulting data 
for C is smooth. Error bars are obtained by determining 
the specific heat from the corresponding finite difference 
as in Eq. (|^) for each sample separately, and computing 
the standard deviation. The error bar is, as usual, the 
standard deviation divided by ^TVsamp — 1- 

In Fig. H the bond energy per spin Ej for two rep- 
resentative L — 8 systems is shown as a function of h. 
For very small values of h all spins point into the same 
direction and so Ej — —3. For large h the spins fol- 
low the random fields and so Ej ^ in this limit. The 
curves in Fig. ^ are stepwise constant functions because 
generically it is not favorable to flip spins if the random 
field is increased by a small amount. However, at certain 
discrete field values, the total energy of another state, 
which differs in the orientation of a cluster of spins, will 




FIG. 3. The upper figure shows the average bond-energy 
[Ej]h per spin as a function a function of the random-field 
strength h for L = 16. The lower figure displays the resulting 
"specific heat", calculated from Eq. (pi) of the text. 



become degenerate with the energy of the ground state 
and for slightly larger values of h the state with the clus- 
ter flipped will become the new ground state. Although 
the total energy is continuous at the field values where 
the ground state configuration changes, the bond energy, 
which is just the first term in Eq. (0), changes discon- 
tinuously. At larger field values the jumps in Ej occur 
closer together and would be difficult to distinguish on 
the scale of a figure. This is why we show, in Fig. ^, data 
for a rather small size. Even for small sizes, the jumps 
occur at different values of h for different samples, and 
so the average value of Ej is expected to be smooth. 

This is illustrated in the upper part of Fig. |^ for i = 16 
which shows a smooth variation of [Ej]h with h. The 
data in the lower part of the figure is the average specific 
heat, obtained as the numerical derivative of the data for 
[Ej]h according to Eq. (||). The specific heat is seen to 
have a peak, as expected. We will investigate the size 
dependence of this peak in Sec. IV. 

In addition to the specific heat, we also calculate the 
susceptibility by considering the response to a small uni- 
form external field H, i.e. we consider the Hamiltonian 



H = — J SiSj — hiSi — H Si . 



(9) 



For each realization, the sign of H is chosen in the di- 
rection of the magnetization of the ground state. This 
prevents the whole system from flipping when applying a 
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FIG. 4. The average magnetization m as a function 
of a uniform external field H near the transition for 
L = A,h = 3.75 (inset: L = 16, /i = 2.8). The solid lines 
represent the results of fits to a parabola, while the dashed 
lines display the tangents at iif = 0; i.e. their slope gives the 
susceptibility. 



magnetic field to a system which is almost ferromagneti- 
cally ordered. The scaling behavior of the magnetization 
should not be affected by this choice. In Fig. |[ the re- 
sult is shown for system sizes L = 4 and L — 16 near 
the values of the field, where the susceptibihty attains a 
maximum. Near H = 0, the data points can be fitted 
very well with a parabola, the coefficient of the linear 
term gives the zero field susceptibility x = dm/dH\H=Q- 
Thus, in order to calculate the susceptibilities, we per- 
form ground state calculations for three different values 
of the uniform field ff„ = nHj^ {n = 0, 1, 2), where, for 
each size, the value of Hl used is shown in Table. |, 
along with the number of samples. We chose the values 
of Hl for each size as follows. For the smaller sizes we 
performed several fields values, as shown in Fig. |4[ to 
determine for what range of fields a parabola accurately 
fitted the data. For larger sizes, finite-size scaling tells 
us that, near the critical point, the characteristic field 
scales with L as L~^" where the "magnetic exponent" 
yH is given by {'J + P)/ v-, with 7 the susceptibility expo- 
nent, and (3 the order parameter exponent. As discussed 
further in Sec. several calculations give /3 ~ 0,7 ~ 2, 
and ;/ ~ 1.3, and so j/jj ~ 1.5. We therefore scale Hl for 
the larger sizes by a factor of roughly L~^-^. 

For each system size, we fit a parabola through the 
three data points for the average magnetization m{Hn)- 
To estimate the error, we performed a jackknife analy- 
sis in which we divided the results for the magne- 
tizations (for each system size and each strength of the 
disorder) into K blocks, calculated the average values K 
times, each time omitting one of the blocks, and then 
performing K fits. The error bar is estimated from the 
variance of the K results for the linear fitting parameter. 
We used if = 50 and checked that the result does not 
depend much on the choice of K . 



L 


^ samp 


Hl 


4 


10^ 


0.05 


6 


60000 


0.025 


8 


40000 


0.016 


12 


30000 


0.008 


16 


23000 


0.005 


24 


27000 


0.0028 


32 


15000 


0.0018 


48 


15000 


9 X 10"** 


64 


9000 


6 X 10"* 


96 


3800 


3 X 10~* 


TABLE I. The maximum number of samples 


A''samp used. 



and sizes of smallest non-zero uniform field Hl, for each sys- 
tem size L. As discussed in the text, the number of samples 
used was larger in the vicinity of the peaks in the susceptibil- 
ity and specific heat than elsewhere. 
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IV. RESULTS 

We have studied random-field systems with sizes from 
L = 4 to L = 96. For each size, simulations were made 
for several different values of h, always averaged over 
many realizations of the disorder. Near the ferromagnet- 
paramagnet phase transition, the number of samples used 
is the largest, ranging from 10^ for the smaller system 
sizes to 3800 for L = 96 for each value h, as shown in 
Table |. With current algorithms, it is in principle pos- 
sible to study even larger system sizes, such as L = 128 
or even L = 256, but, using the LEDA algorithms, these 
need more memory than the 512 MBytes available to us. 
Hence we have restricted our study to L < 96, which is 
still much larger than sizes that can be simulated using 
Monte Carlo simulations. 

In the thermodynamic limit the singular part of the 
specific heat diverges according to 



Cs K A±\h~ he 



(10) 



where the amplitudes and j4_ refer to h > he and h < 
he respectively, and a is the specific heat exponent. In 
addition there is a regular piece of the specific heat, Ciog, 
which is finite at the critical point and so dominates there 
if a < 0. In a finite system, finite-size scaling predicts 
that 



(11) 



where v is the correlation length exponent. The specific 
heat peak will occur when the argument of the scaling 
function C takes some value, oi say, so the peak position 
h*{L) varies as 



h*{L) — he ~ aiL 



(12) 



and the value of the singular part of the specific heat at 
the peak varies as 



ol/u 



(13) 



In Fig. H the specific heat C is shown as a function 
of the random-field strength h for selected system sizes. 
The error bars are obtained from the standard deviation 
of the data for different samples, and are quite small be- 
cause a large number of samples have been averaged over, 
see Table |. A clear peak can be seen, which moves to 
the left and increases in height with increasing system 
size. The number of samples used is larger near the peak 
to compensate for the greater sample to sample fluctua- 
tions in this region. For each system size, we performed 
parabolic fits to the region of the peak to obtain h*{L) 
and the height of the peak, C""^'^(L). The shift of the 
maximum according to Eq. ( p^ can be used to estimate 
the infinite-size critical strength of the random field, he 
and the correlation-length exponent v. The best fit gives 




FIG. 5. "Specific heat" C, calculated from Eq. (g), as 
a function of the random-field strength h for system sizes 
L = 4, 8, 16, 32, 64 and 96. The vertical dashed line indicates 
the location of the critical value of the random field, he = 2.28, 
see Eq. . The inset is an enlargement of the peaks for the 
larger sizes. 



see Fig. O. We determined the probability Q that the 



value of 



EL( 



'yi-f(xi)\2 



with N data points 



{xi,yi ± (Ti) fitted to the function /, is worse than in 
the current fit to quantify the quality of the fit. Here 
we get Q = 0.20, which is fair. 

Next we try to estimate the specific heat exponent by 
looking at how the peak value C"^'^^ scales with L. If a = 
one expects logarithmic divergence and the simplest 
hypothesis is to fit the data to 



blogL, 



(15) 



where the constant term a comes partly from the regular 
piece of the specific heat. However, Fig. ^ shows that 
this does not work. A plot of C"^'^^ against L (on a log 
scale) shows clear curvature, suggesting that the height 
of the specific heat will saturate to a finite value as L 
increases. If one considers only the data points for sizes 
L = 4, . . . , 16, as in Ref. 1 13|, a negative curvature is still 
visible, but the result is much less clear. 

A peak height which saturates for L — > oo implies that 
a is negative, in which case the specific heat has a finite 
cusp at the critical point, rather than a divergence. We 
have therefore tried a fit of the form 



/ic = 2.28 ± 0.01, l/i/ = 0.73 ± 0.02, 



(14) 



(16) 
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FIG. 6. A plot of the random field where the specific heat 
attains its maximum, as a function of system size L. The 
solid line shows a fit to the function h*{L) — he + aiL"^^" 
with he = 2.28, 1/v = 0.73, and a\ = 2.55. The inset shows 
the data as a function of L~^^'^ . 

in which Coo comes from the regular part of the specific 
heat, yielding, 

Coo = 2.84 ±0.05, a/z^ = -0.48 ±0.03. (17) 

This fit is shown in the inset of Fig. |^. The quality 
of the fit, Q = 0.05, is not very good. We have tried 
different fits using only the larger system sizes, which 
increases the quality of the fit slightly, but the resulting 
error bars are very large. The central estimate for a 
actually becomes more negative if we only include the 
larger sizes. The rather poor fit may indicate difficulty in 
accurately estimating the error bars for the location and 
height of the specific heat peak. Our analysis suggests 
that the specific heat exponent is strongly negative, in 
agreement with Rieger |I3| though we cannot rule out a 
leading singularity with a ~ and a sufficiently small 
amplitude that it is hard to see in our data. 

To look for this possibility, we also tried more compli- 
cated fits including corrections to scaling of the form 

= Coo + a2i"/"(l + bL-'^), (18) 

where m is the leading correction to scaling exponent. 
The data did not determine all the parameters cleanly, 
and the fit program which works iteratively, con- 
verged to different results depending on the starting val- 
ues, and whether any of the parameters were held fixed. 



The solutions we found were of two types: (i) the fit is 
the same as that in the simpler fit of Eq. ( [iq ) (i.e. w is 
essentially zero and a/v and the other parameters are the 
the same as found in the simpler fit), (ii) lj is quite small, 
02 is very large, and b is negative such that 1 + bL^^ is 
close to zero. Thus, in the second type of fit, the data 
is represented as two singularities with large amplitudes 
which almost cancel. This does not seem physical. The 
fitting routine did not converge to a solution with a lead- 
ing singularity which has a small amplitude and a ~ 0, 
plus a correction term with a larger amplitude. 

We will discuss our specific heat results further in 
Sec.0 
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L 

FIG. 7. The maximum C'^'^^ of the specific heat as a func- 
tion of system size L with logarithmically scaled L-axis. The 
dashed line is a tangent to the data and a comparison be- 
tween it and the data demonstrates that (7™^" grows slower 
than logarithmically with system size. The solid line shows a 
fit to the function C""^'=(L) = Coo + a2L°'''' with Coo = 2.84, 
a/u = —0.48 and a2 — —3.52. The inset shows the data and 
the fit as a function of L"''". 

The susceptibility x a function of h is presented 
in Fig. H for selected system sizes. It is seen that the 
height of the peak grows much faster than for the specific 
heat. To analyze the divergence of we have again fitted 
parabolas to the data points near the peak to obtain 
the positions h*{L) and x^^^{L) of the maximum. By 
fitting the data for L > 32 to a function x™'"'(i) = 
a^L^^^, where ry describes the decay of the "connected" 
correlations at criticality, we obtain (Q = 0.63) 

= 0.50 ±0.03, (19) 
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FIG. 8. Susceptibility x as a function of the random-field 
strength h for system sizes L — 8, 16, 32, 64, and 96. Only 
data near the peaks is shown because the data away from the 
peaks had lower precision. 
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FIG. 9. The maximum x™^" of the susceptibility as a func- 
tion of system size L in a double logarithmic plot. The solid 
line represents a fit to the function x^'^^iL) = a3L^~^ , for 
sizes L > 32 yielding 2 — rj — 1.50 and 03 = 0.095. 



see Fig. |[ 

Finally, we have also estimated he and the correlation- 
length exponent from the susceptibility data using Eq. 
(p^), as we did for the specific heat. Using only sizes 
L > 32 (Q = 0.84), we find 



/ic = 2.29 ±0.01, l/zy = 0.81 ±0.05. 



(20) 



jrees with that obtained from the 
) , while the estimate for l/v dif- 
h by slightly more than the sum 
of the error bars, probably indicating some systematic 
corrections to scaling. 



This estimate of he a£ 
specific heat, see Eq. (uA 
fers from that in Eq. (fiA 



V. DISCUSSION 

We have determined the "specific heat" of the random 
field Ising model at T = using optimization algorithms. 
The height of the peak increases less fast than logarith- 
mically with system size, and a finite-size scaling analysis 
gives the exponents shown in Eqs. and (p^). From 
the analysis of the susceptibility, the exponents shown 
in Eqs. ( p^ ) and ( pO| ) are obtained. The final results we 
quote are 



hc^ 2.28 ±0.01, z/ = 1.32 ±0.07 
a =-0.63 ±0.07, 77 = 0.50 ± 0.03. 



(21) 



To determine v and its error we have taken both the val- 
ues in Eqs. ( p^ ) and ( |20| ) and used the difference between 
them as a measure of the systematic error. The errors 
for 77 and he are purely statistical. The error for a comes 
both from the error in 1/ and the statistical error in a/u. 

Our results for he are compatible with the values 2.29± 
0.04 1^], 2.26 ±0.01 m, and 2.270 ±0.005 ]l| obtained 
from ground-state calculations of systems of similar size. 
Values of v obtained from ground-state calculations are 
1.37±0.09 |12) and 1.19±0.08, (2|], which agree well with 
our result. Ref. |2^ argued for a first-order transition, 
but assuming scaling with respect to the field, a value of 
h' = 1.25 ± 0.06 was estimated, also in agreement with 
our result. However, if a power law correction to scaling 
was taken into account, instead the result 1.52 (without 
error bars) was found. 

The scaling exponent 77 describing the susceptibility, 
has not been obtained from exact ground-state calcula- 
tions so far. In a Monte-Carlo simulation ||l^ a value 
of 0.50 ± 0.05 was found, which is compatible with our 
result. 

The most significant result of this paper is that for the 
specific heat, namely a — —0.63(7). This agrees well 
with the values a/v = -0.45 ± 0.05, j/ = 1.1 ± 0.2 found 
by Ref. || and a = -0.55 ± 0.20 found by Ref. H, 
both using Monte Carlo simulations on small systems. 
However, as we shall now see, it appears inconsistent with 
values for other exponents and expected scaling relations. 
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At conventional second-order phase transitions, all ex- 
ponents can be related to two (e.g. v and 77) by scaling 
relations. However, because the fixed point of the RFIM 
is at T = with temperature a "dangerous irrelevant 
variable", a modified set of scaling relations has been 
proposed [^j^Js^,^ , which involve three independent ex- 
ponents. Scaling relations which do not involve the space 
dimension, e.g. 



(22) 



are unchanged, but "hyperscaling" relations involving the 
space dimension d, have d replaced by d—9, where 6, the 
third exponent, is the scaling exponent for the temper- 
ature at the fixed point. An example of a hyperscaling 
relation which is relevant to the specific heat is 



{d-e)v = 2-a. 



(23) 



Gofman et al. have proposed that the Schwartz-Soffer 
PH inequality, which can be expressed as 77 > 2 — 0, is an 
equality, in which case there are only two independent 
exponents again (though the hyperscaling relations are 
different from those in conventional two-exponent scal- 
ing). Our results are consistent with this, since /3 ~ 
implies that 9 ~ 1.5, see e.g. Ref. ||l2[, and we have 
already found that is about 0.50, see Eq. (|l9|). 

26 1 (3 (the most accu- 
[ ) , and our value for 



Other works have found 
rate value is 0.017 ± 0.005 in Rcf. ^ 
7, obtained from 7 = (2 — 77)1/ is about 2.0 in agreement 



with series expansion work of Gofman et al. |40|. Hence 
Eq. (^) predicts a ~ 0, quite different from the value of 
about —0.63 that we find by direct calculation. 

As noted above, the result /3 ~ implies that 6 ~ 1.5, 
so Eq. ( p^ ) gives a ~ 2 — l.^v. Using our value oi v — 
1.32 ± 0.07 this yields a = 0.0 ± 0.15. In other words, 
Eq. ( p3| ) also predicts that a is close to zero. 

We have seen that the two scaling relations above 
would be consistent if we inserted a ~ 0, which is the 
experimental value However, by direct calculation, 
we obtain a strongly negative result, a ~ —0.63, consis- 
tent with earlier work [|l3| on much smaller sizes. Thus 
the problem with the value of the specific-heat exponent 
has now been strongly reinforced by our calculations on 
much larger lattices. 

Possible explanations for this discrepancy are: 

• The specific heat diverges but slower than logarith- 
mically. Examples of this, which are known to oc- 
cur in other systems, are a fractional power of a 
log and a log-log variation. However, there are no 
calculations which predict this type of behavior for 
the RFIM. Furthermore, attempts to fit our data 
to this type of behavior were not very successful. A 
related possibility, which does not seem impossible 
looking at Fig. ^, is that a = might be realized 
by a jump in the specific heat, with a lower value 
in the ferromagnetic region, the opposite of what 
occurs in mean field theory. 



• The regular contribution to the specific heat varies 
rapidly near the critical point. Since /3 ~ 
the magnetization increases very rapidly below he 
(leading to the very rapid drop in the specific heat 
seen in Fig. ||) . If much of this drop comes from the 
regular part of the specific heat it would be difficult 
to extract the singular part. 

• There are very strong singular corrections to finite- 
size scaling which leads to the most singular term in 
the specific heat being numerically small compared 
with correction terms, even for the quite large range 
of sizes that we have studied here. If there are 
strong corrections to scaling, perhaps the values of 
other exponents, in addition to a, could be affected 
too. 

• Scaling does not hold. We find this possibility to 
be the least palatable. 

Since /? ~ 0, it is interesting to ask whether the tran- 
sition might be first order and whether this might be the 
origin of the surprising value of a. The transition at low- 
T is first order in mean field theory for field distributions 
with a minimum at zero field |4^ . A first order transi- 
tion for Gaussian distribution has also been suggested 
for dimension less than four based on series expansion 
work fist . If the transition is first order, it must be very 
weakly so, since fluctuation effects are very large. Fur- 
thermore, one would then expect a latent heat, which, 
in a finite-size system, gives a specific heat diverging as 
the volume L'^. In our results, we do not see any diver- 
gence, let alone a strong one like this. In addition, the 
most detailed numerical study jl^] claims that /3 while 
very small, is greater than zero. Even if the transition 
were ultimately first order, the effective exponents found 
should be those of the close-by second order transition, 
and so should satisfy scaling. We therefore don't feel that 
the possibility of a first order transition explains why our 
value for a does not satisfy scaling. 

In addition to critical exponents, it is useful to dis- 
cuss amplitude ratios, since these are also universal, see 
Ref. iQ and references therein. For the specific heat 
amplitudes, and A^, defined in Eq. (|l0|), one can 
show pa that A^/A^ = 1 for a logarithmic divergence 
(a = 0). Furthermore, for n-component models with- 
out random fields one has A^ /A- > 1 for a < 
and A+/A- < 1 if a > 0. This implies that, for both 
signs of a, the specific heat decreases from its peak faster 
on the paramagnetic side than on the ferromagnetic side 
(we are grateful to D. Belanger for pointing this out). 
By contrast, the situation is reversed in our data, see 
Fig. H where the specific heat appears to decrease faster 
for h < he- Whether this indicates that the amplitude 
ratio is very different in the presence of random fields, or 
that corrections to scaling are large compared with the 
leading singularity for this range of sizes remains to be 
seen. 
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Clearly more work is needed to understand the spe- 
cific heat of the RFIM. Since several recent large-scale 
numerical calculations, including ours, have used fairly 
sophisticated algorithms, it is unlikely that a numerical 
breakthrough is imminent. Hence a better theoretical 
understanding, especially of corrections to scaling, will 
be needed to sort out this problem. 

Additional Note: After this work was submitted we 
received the final version |Q of Ref. in which, mo- 
tivated by our work, they computed the bond energy 
using ground state methods. They did not numerically 
differentiate the data to get the specific heat but directly 
analyzed data for the bond energy at the bulk critical 
field, the dashed line in Fig. |^. The size dependence in- 
volves the exponent (1— a) /i' from which they find results 
compatible with a — 0. That they get a different result 
from ours by, in effect, considering a different region of 
the scaling function in Eq. (|ll|), indicates that there are 
large corrections to finite size scaling even for such large 
sizes, or possibly that a ~ corresponds to a discon- 
tiniuty in the specific heat. Both these possibilities were 
discussed above. Further work is needed to clarify the 
situation. 
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